Spinodal fractionation in a polydisperse square well fluid 
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Using Kinetic Monte Carlo simulation, we model gas-liquid spinodal decomposition in a size- 
polydisperse square well fluid, representing a 'near-monodisperse' colloidal dispersion. We find that 
fractionation (demixing) of particle sizes between the phases begins asserting itself shortly after the 
onset of phase ordering. Strikingly, the direction of size fractionation can be reversed by a seemingly 
trivial choice between two inter-particle potentials which, in the monodisperse case, are identical - 
we rationalise this in terms of a perturbative, equilibrium theory of polydispersity. Furthermore, 
our quantitative results show that Kinetic Monte Carlo simulation can provide detailed insight into 
the role of fractionation in real colloidal systems. 

PACS numbers: 82.70.-y, 64.75. Gh, 64.60.-i 
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I. INTRODUCTION 

Colloidal dispersions consist of particles of mesoscopic 
size dispersed in a fluid solvent, a definition that em- 
braces a huge variety of natural and synthetic substances. 
Milk, blood, viruses, paints, cosmetics and even radioac- 
tive waste are all colloidal, in that they are made up 
of some continuous background solvent in which one or 
many species of colloidal particle are dispersed. Since the 
microscopic dynamics of the colloidal particles are dom- 
inated by thermal fluctuations in the solvent, they can 
be studied using familiar thermodynamic and statistical 
mechanical approaches Therefore, as well as being 
studied for their own sake, colloids are seen as a use- 
ful analogue of ordinary molecular fluids and - especially 
due to the relatively accessible length and time scales in- 
volved [2| - as an experimental testing ground for more 
general theories 0,0] • Nevertheless, the analogy between 
molecular and colloidal fluids is not exact. One key differ- 
ence is that colloidal particles are inevitably polydisperse, 
i.e. they exhibit continuous variation in particle proper- 
ties such as size, charge or shape. Polydispersity leads to 
complex phase behaviour and complicates the theoreti- 
cal and experimental treatment of complex fluids [sT-flci] . 
Increasing the breadth, precision and utility of our un- 
derstanding of soft matter physics thus relies in part on 
a deeper understanding of the behaviour of polydisperse 
systems. 

Much existing work has concerned itself with the equi- 
librium thermodynamic effects of size-polydispersity pi- 
ll, ITlTfl6l ] . Polydispersity is predicted and/or observed 
to cause destabilisation of the solid phase, fractionation 
(demixing) of particle sizes between phases [lj]> EH and 
an apparent 'terminal polydispersity' beyond which the 
system must split into multiple solid phases to reach equi- 
librium Q. 

In contrast, theoretical work on the dynamics of poly- 
disperse phase behaviour [l?], EH is rare, and the dynam- 
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ical picture of how polydisperse systems actually do or 
do not arrive at equilibrium is far from complete [19l - l23l | . 
Moreover, in such studies as do exist, polydispersity is 
often introduced only to inhibit crystallisastion so that 
high density amorphous states can be studied [13, HH ■ 

Fractionation is thought to be involved in the nucle- 
ation of crystals from polydisperse fluids (2l|-|23|. with 
frustration of this mechanism held responsible for the 
experimental inaccessibility of some complex equilibrium 
phase diagrams Q. In practical applications, fractiona- 
tion of particle species between separating phases may be 
intended, for example to purify a mixture, or unintended. 
In either case, it is important to consider consequent ef- 
fects on related properties (composition, crystal lattice 
parameter etc.) in the resultant phases. 

In the present work we address the issue by dynam- 
ically simulating gas-liquid spinodal decomposition in a 
size-polydisperse colloidal fluid. We observe fractiona- 
tion phenomena on simulation timescales, and uncover a 
surprising example of the sensitivity of such phenomena 
to the details of particle interactions, which we explain 
by appealing to equilibrium theory. 



II. SIMULATION SETUP 

The system studied in this work consists of either 
N = 30000 or N = 4000 spherical particles whose diam- 
eters d n are drawn from a Bates (pseudo-Gaussian) dis- 
tribution of polydispersity (defined as ratio of the stan- 
dard deviation to the mean) cr = 0.06 or er = 0.2 [26| . 
The simulation cell is cubic with periodic boundaries. 
We set the length unit as the mean hard core diameter 
(d) . In order to mimic the qualitative features of any at- 
tractive interaction in the colloid, the hard particle cores 
are surrounded by a pairwise additive square well po- 
tential of (mean) range A = 1.15 and depth u = 1.82, 
where the energy unit is ksT. In the monodisperse limit, 
these parameters and the overall 'parent' volume frac- 
tion 4> p = 0.229 place the s yste m in the middle of the 
gas- liquid coexistence region 27| , where phase separation 
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would be expected to proceed by spinodal decomposition. 
For the value of A used here, the resulting gas-liquid co- 
existence is metastable with respect to fluid-solid phase 
separation which, on the timescale simulated, we do not 
observe. 

Note that there is a choice to be made in precisely how 
the specific attraction range between any two particles i 
and j should depend upon their size. We have simulated 
two possibilities. The 'scalable' dependence is motivated 
by its simplifying effect on the relevant theory viz. the 
equilibrium effects of fractionation [l2| : 



The Metropolis acceptance criterion, popular for equi- 
librium Monte Carlo (MC) simulations, assigns an equal 
probability (unity) to all moves which do not increase the 
system's energy, even those which decrease it. In order 
to more naturally represent the dynamics of attractive 
particles, we instead employ an acceptance probability 



exp(-£ new ) 



accept 



exp(-E new ) + exp(-.E old ) 



(3) 



^scal(r) 



oo if r < dij 
—u if d^ < r < Xdij 
if r > \da 



(1) 



where d^ = {d. L +dj)/2. The 'non-scalable' dependence 
is chosen for its similarity to depletion attractions [28j j . 
in which the attraction range depends on the hard core 
diameter via an additive constant, e.g. the radius of gy- 
ration of a dissolved polymer: 



oo if r < dij 

Kon-scalM = { ~U if dij < T < d tj + (A - 1) 

if r > d^ + (A - 1) 



(2) 



In the monodisperse case, dij = 1 for all particle pairs so 
that the two definitions are strictly identical. Even in the 
polydisperse case the mean attraction range is unaffected 
by the choice of dependence, but it will nevertheless turn 
out that this easily-overlooked detail can make a radical 
difference to the fractionation of the system. 

Kinetic (or 'dynamic') Monte Carlo (MC) simulation 
is distinct from equilibrium MC methods in that it is re- 
stricted to dynamically realistic trial moves, i.e. small 
stochastic 'hops' representing the movement of particles 
performing Brownian motion. This is in contrast to equi- 
librium algorithms in which unphysical moves (e.g. clus- 
ter rearrangements or biased sampling) may be used to 
speed up equilibration so long as detailed balance is re- 
spected. For small step sizes, dynamical results from 
Kinetic MC simulation compare well with those of Brow- 
nian dynamics [29| . The simplicity of Kinetic MC and its 
natural suitability for representing inertia-free stochastic 
motion has resulted in a number of applications in soft 
matter simulation [24|, [3(1 HH . 

In our implementation, trial particle moves are drawn 
from a Bates distribution centred on zero, and the mean 
step size for particle i is Ay/(d)/di, so as to give the 
correct dependence of the short time Stokes-Einstcin dif- 
fusion coefficient upon particle diameter. (The step pa- 
rameter is set to A = 0.02, and reducing this value does 
not alter the results presented here). The time unit is 
the time td taken for a free particle of diameter (d) to 
diffuse a distance equal to its own diameter. 



to move from a configuration with energy E Q \^ to one 
with E ncw . Moves that conserve energy are now accepted 
with probability 0.5, so that moves that bring a particle 
into the reach of attractive wells can be assigned a higher 
probability, representing more rapid movement down a 
potential gradient. 

The Kinetic MC algorithm allows large systems to be 
simulated on a physically relevant timescale, necessary 
for measurement of the required statistics for fractiona- 
tion. However, as in comparable Molecular or Brownian 
dynamics algorithms, hydrodynamic interactions (HI) 
are neglected. Given the importance of HI in spinodal 
decomposition, this requires some comment. In an ap- 
proximate theoretical scheme, for example, HI have been 
shown 34] to enhance longer-wavelength decomposition 
while suppressing shorter wavelengths. Measuring the ef- 
fect of such corrections on fractionation in a polydisperse 
system would require for instance explicit treatment of 
the solvent via Molecular dynamics, or a combined Brow- 
nian dynamics/lattice Boltzmann method, considerably 
increasing computational expense. 

Nevertheless, we expect that the phase separation ob- 
served here is at least qualitatively representative of spin- 
odal decomposition in a real system and, particularly 
given that HI can be screened (e.g. by direct interac- 
tions in dense systems [35]), constitutes a relevant and 
important baseline in the low-HI limit. Given this, and 
since our objective is to collect good statistics for the ob- 
servation of fractionation, we leave the incorporation of 
HI to future work. 



III. SPINODAL DECOMPOSITION 

The system is initialised as an amorphous fluid, where- 
after the square well attraction is turned on (defining an 
instantaneous quench at t — 0) and the system evolved 
according to the Kinetic Monte Carlo algorithm outlined 
above. In addition to direct observation through time 
(FIG. [1]), the development of a low-g peak in the struc- 
ture factor S (q) is used to track the gas-liquid phase sep- 
aration (FIG. [5]). As shown in FIG. O this peak has 
already started to shift to lower q after 2 or 3 tj, indi- 
cating coarsening of the gas-liquid texture. To confirm 
the spinodal nature of the separation, FIG. [3] shows the 
integrated area under the peak increasing on a timescale 
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FIG. 1. Illustrative simulation snapshots at (clockwise from 
top- left): t = 0, t = 4, t — 64 and t — 760 for a slice 3 units 
wide through the system, showing the coarsening spinodal 
texture. 



of hundredths of td- We detect no 'lag' period associated 
with nucleation-driven phase separation, so we conclude 
that the phase separation is indeed spinodal in nature, 
as expected for our parameters 
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Other dynamical observations, to be presented else- 
where, are in agreement with Ref. [36] in suggesting 
that, as in the monodisperse case, the gas-liquid coex- 
istence resulting from the spinodal decomposition is in 
turn metastable with respect to the growth of a crys- 
talline phase that would coexist with a tenuous vapour. 
On the timescales observed here and in the absence of a 
template to seed nucleation, we do not observe crystals. 




FIG. 2. Structure factor S(q) for one simulation trajectory, 
averaged over a) t — 1 — 2 and b) t = 3 — 4. The usual 
fluid peak is present at q ~ 7, while a peak corresponding 
to the gas-liquid texture appears at q m 2 and shifts to a 
lower wavenumber as it grows, indicating coarsening of the 
two phases. 



IV. FRACTIONATION 

To check for fractionation, a per-particle distinction 
between the gas and liquid phases is required. We count 
a particle as liquid if it has more than 12 neighbours 
within 2 length units, this threshold corresponding to a 
local number density in the middle of the gas-liquid bin- 
odals. Throughout the simulation the mean diameters of 
particles in the gas and liquid regions, (d) gas and (d)u qu id) 
are recorded. 

We present fractionation data from 4 large simula- 
tions with N — 30000, where 2 independent initial con- 
figurations have been simulated according to both the 
scalable and non-scalable definitions of the inter-particle 
potential. In addition, we ran 20 smaller simulations 
(N = 4000), 10 for each potential, to confirm the qualita- 
tive nature of the fractionation effects, check for sensitiv- 
ity to initial configuration or system size, and investigate 
longer timescales. Finally, having confirmed fractiona- 
tion in the near-monodisperse regime, we ran 2 simula- 



tions at a higher polydispersity a — 0.2 and N = 4000 
in which we expect the strength of fractionation, scal- 
ing with a 2 [12j . to increase by an order of magnitude 
relative to the other simulations. 

FIG. |4] shows the evolution of (d) gas and (rf)u qu id in 
the large system for the scalable and non-scalable po- 
tentials V sca \(r) and V uon - sca x(r) . Some fractionation of 
particle sizes is apparent during the first few hundred 
time units. As shown in FIG. [5j the fractionation in the 
smaller, longer-time simulations agrees within error with 
the data from the larger simulations. Thereafter, the 
fractionation in these simulations shows no detectable 
change beyond t « 2500 up to the maximum time simu- 
lated. FIG. IB] confirms the same qualitative fractionation 
at higher polydispersity, with the difference in mean size 
between the phases increased by approximately an order 
of magnitude as compared to the a = 0.06 case. 

Since fractionation relies on self-diffusion of individual 
particles between the phases, whereas bulk phase sepa- 
ration takes place via faster, collective rearrangement of 
particle density [l7J , it might be expected that the speed 
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FIG. 3. The area under the Iow-g peak, found by integrat- 
ing from q — 0.26 to q = 2.56, averaged over 3 indepen- 
dent quenches, for times immediately after the quench. Al- 
though background noise, indicated by the value at t = 0, 
is strong, the area under the peak is clearly increasing from 
these very early times, consistent with spinodal decomposi- 
tion. The dashed line is a guide to the eye. 



of spinodal decomposition would negate any possibility 
of fractionation in the early stages of phase separation, 
if the particles could not rearrange to the required de- 
gree faster than the spinodal texture coarsens. Therefore, 
while our results leave open the question of further 'ripen- 
ing' of the phase composition on much longer timescales, 
it is significant that fractionation is apparent almost im- 
mediately after spinodal decomposition begins. 

More strikingly, there is a qualitative dependence on 
the choice of inter-particle potential - the non-scalable 
potential causes the liquid to incorporate smaller par- 
ticles, whereas the scalable potential causes the smaller 
particles to be incorporated into the gas. Recall that 
Kicai(0 and V^n-scaiM nave the same well depth and 
mean range - they differ only in the precise dependence 
of their range upon the mean diameter dij of the particles 
between which the potential is acting. The two poten- 
tials coincide by definition in the monodisperse case (for 
which they produce exactly identical simulation trajec- 
tories), but are 'split' by a mild degree of polydispersity 
so as to produce opposite fractionation effects. 

We can explain this surprising result with an equilib- 
rium theory of polydispersity [12[ ■ While the system be- 
ing simulated here is clearly not at equilibrium, it is rea- 
sonable to invoke such a theory in order to give some in- 
dication of what metastable state the system is relaxing 
towards. In Ref. [HJ, perturbation of a monodisperse ref- 
erence system is used to predict equilibrium properties of 
polydisperse systems, requiring as input only properties 
of that monodisperse reference. The approach becomes 
exact in the limit of weak polydispersity, i.e. if (e 2 ) is 
small where e is a particle's deviation from the mean size 
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FIG. 4. (colour online) (d) gas (black) and (d)ii qu id (red) 
through time for the two definitions of the inter-particle po- 
tential, for 4 simulations with N = 30000. The parent mean 
particle diameter is (d) = 1. a) For the scalable potential 
K, C ai(r) the liquid prefers on average larger particles than the 
gas. b) The non-scalable potential V uou . a ca.i(r) causes the liq- 
uid to prefer smaller particles. Error bars represent the stan- 
dard error on the mean. 



(d), in units of (d). The applicability of that theory is 
the reason for our primarily focusing on the low-er regime 
in the simulations. 

For fractionation, the key thermodynamic potential is 
A(p) = pdpf x (e)/de, a function of number density p de- 
scribing how the excess chemical potential for a particle 
in a monodisperse system varies with e. Once A{p) and 
the polydispersity a of the parent distribution are known, 
the fractional difference in mean particle size between the 
phases is given by: 



-[A/p\y 



(4) 



where the labels I, g indicate quantities evaluated at the 
metastable coexistence densities of the liquid and gas 
phases respectively. Thus, if the difference in Aj p be- 
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FIG. 5. (colour online) As FIG. [4] for 20 simulations with 
N — 4000, run for longer times. The fractionation shows no 
significant change beyond t ~ 2500, and is consistent in its 
magnitude and direction with the data in FIG. [4] 



tween the two phases, [A/p] 1 , is negative, the predicted 
mean size in the metastable liquid phase will be greater 
than that in the gas. The a 1 dependence is consistent 
with the difference of around one order of magnitude in 
the fractionation between the a = 0.2 and a = 0.06 sim- 
ulations. 

To find A, an approximate free energy density / is re- 
quired for the monodispcrse square well fluid. We use the 
attractive part of the interaction potential U(r) = —u as 
a perturbation to a hard-sphere reference system with the 
Percus-Yevick radial distribution function t/Hs( r )i giving 



/hs + y 



47rr 2 gns( r )V (r) dr 



(5) 



where /hs is the Carnahan-Starling hard-sphere free en- 
ergy density [13] and the integral is over the square well 
range. 

For V Bca i(r), finding [A/p] 1 is made simpler because 
the definition is such that all lengths, including the range, 
scale with the hard core radius. In this special case, it can 
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FIG. 6. (colour online) As FIG. [H for higher polydispersity 
(<j = 0.2) simulations. The direction of fractionation is the 
same as in the lower-<r simulations, and the difference in mean 
size between the phases is increased by approximately an or- 
der of magnitude. 



be shown [l2[ that A = 3(P — p) where P is the osmotic 
pressure calculated from the free energy density. FIG. [7] 
shows that, for the scalable potential, A/p decreases 
across the gas-liquid coexistence region (</> fl w 0.05 to 
4>i w 0.40) so that [A/p] 1 < 0, correctly predicting the 
observed fractionation of larger particles into the liquid 
phase in FIG.|Ua). For Vnon-scaiM, A/p increases mono- 



tonically, [A/p] l g > and the liquid contains on average 
smaller particles, in qualitative agreement with FIG. [4] 
b). 

This is intuitively reasonable: from Equations [T] and 
[2l it is clear that the range of V sca \(r) increases more 
strongly with dij than that of K ao n-scal('")- I n the non- 
scalable case, increasing the diameter of particles in the 
liquid results in a greater increase in chemical potential 
than in the gas - the incorporation of larger particles 
into the liquid is unfavourable. In the scalable case, the 
marginally stronger dependence of the well range on d^ is 
such that an increase in particle diameter reduces chem- 
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FIG. 7. A/p shown as a function of volume fraction cj> for 
a monodisperse reference system with u = 1.82 (solid lines), 
used to predict the fractionation of particle size at gas-liquid 
coexistence Dashed lines show u = 1.74 for comparison. 
The lower two curves are for the scalable interaction potential 
V^cai(r), while the upper two curves correspond to the non- 
scalable potential V non . sca i{r) . 



ical potential, as more particles can fit into the square 
wells - this effect is enhanced in the liquid, compared 
with the gas, due to its smaller inter-particle spacing, so 
that the liquid now favours larger particles. 

For a quantitative comparison, we use Equation U to 
estimate the extent of 'equilibrium' fractionation at gas- 
liquid coexistence by evaluating [A/ p] l g from the curves 
in FIG. [7] Taking the coexistence volume fractions from 
27], (j) g « 0.05 and <j>i 0.40, we find [A scal / ' p] l g = -5.31 
and [A non . sca x/ ' py g = 2.20, giving for the a = 0.06 system 



[{d) sc J g = 0.019, [(d) non . sc J g = -0.0079 . (6) 



It appears that, while the fractionation observed on the 
simulated timescale is in the direction of equilibrium, it 
is substantially weaker than that required to minimise 
the free energy of the metastably coexisting phases, con- 
sistent with the hypothesis that the compositions of the 
phases would typically relax far more slowly than their 
overall densities [17J. Given this, it is all the more re- 
markable that the qualitative sensitivity of fractiona- 
tion to the distinction between Kj Ca i(r) and V non . sca \(r) is 
present even at the very early stages of phase separation. 



V. CONCLUSION 

Whilst dynamical simulations of crystallisation in 
bidisperse mixtures have previously found evidence of 
fractionation [2(| these are, to our knowledge, the first 
such measurements in the early-stage phase separation 
of a truly polydisperse model colloidal fluid. The equi- 
librium theory used here has been tested in experimen- 
tal scenarios [ll|> HU, 39] in which fractionation measure- 
ments are typically taken after a long period of equili- 
bration. For a typical colloidal particle of radius « 100 
nm in water, the diffusion time is td ~ 10~ 3 s in which 
case the simulation time used here represents the first 
few seconds of phase separation. Experimentally, mea- 
surements in these early stages are difficult or impossible 
in the absence of large single-phase regions or when the 
system is evolving quickly. 

The data presented show that the beginnings of frac- 
tionation can take effect from the outset of spinodal de- 
composition, but with a smaller magnitude than pre- 
dicted by equilibrium calculations. Although we fo- 
cused on the low-polydispersity regime in which the per- 
turbative theory used here becomes exact, simulations 
at higher polydispersity confirmed the same qualitative 
physics, accompanied by an order of magnitude increase 
in the difference in mean particle size between phases, 
as expected from Equation 01 We have also found an 
example of the subtle interplay between polydispersity 
and other system parameters, in that the choice between 
two inter-particle potentials which are identical in the 
monodisperse case can qualitatively change the predicted 
and observed direction of fractionation in a slightly poly- 
disperse system. 

The latter finding starkly illustrates the care that must 
be taken as the treatment of polydisperse systems be- 
comes more sophisticated: under a 'near-monodisperse' 
approximation, the choice between the two inter-particle 
potentials studied here might have been expected to have 
a negligible effect on the results and therefore be made 
simply to maximise theoretical expediency. On the con- 
trary, we have shown that this choice can dramatically 
affect (and in this case, reverse) the fractionation ob- 
served, a fact that must be taken into account if precise 
knowledge and control of the composition and dependent 
properties of the two phases is desired. Further, the fact 
that dynamical simulation can resolve such phenomena 
highlights the detailed role it may play in bridging the 
gap between equilibrium theories of colloidal polydisper- 
sity and the behaviour of real, dynamical systems. 
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